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Abstract 

A spectral method is described for solving coupled elliptic problems on an interior and 
an exterior domain. The method is formulated and tested on the two-dimensional interior 
Poisson and exterior Laplace problems, whose solutions and their normal derivatives are 
required to be continuous across the interface. A complete basis of homogeneous solutions 
for the interior and exterior regions, corresponding to all possible Dirichlet boundary val- 
ues at the interface, are calculated in a preprocessing step. This basis is used to construct 
the influence matrix which serves to transform the coupled boundary conditions into con- 
ditions on the interior problem. Chebyshev approximations are used to represent both the 
interior solutions and the boundary values. A standard Chebyshev spectral method is used 
to calculate the interior solutions. The exterior harmonic solutions are calculated as the 
convolution of the free-space Green's function with a surface density; this surface density 
is itself the solution to an integral equation which has an analytic solution when the bound- 
ary values are given as a Chebyshev expansion. Properties of Chebyshev approximations 
insure that the basis of exterior harmonic functions represents the external near-boundary 
solutions uniformly. The method is tested by calculating the electrostatic potential resulting 
from charge distributions in a rectangle. The resulting influence matrix is well-conditioned 
and solutions converge exponentially as the resolution is increased. The generalization of 
this approach to three-dimensional problems is discussed, in particular the magnetohydro- 
dynamic equations in a finite cylindrical domain surrounded by a vacuum. 

Key words: influence matrix, spectral method, Chebyshev polynomials, boundary integral 
method, magnetohydrodynamics. Green's functions, harmonic functions, Laplace's 
equation, exterior problem 



1 Motivation 

The search for a self-sustaining magnetohydrodyamic dynamo has taken on great 
momentum in recent years, as researchers have sought to produce dynamos in the 
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laboratory [I1I2I3I4I5II and in simulations II6I7I8I9I10I11I12I13I15I16II . One of the 

fundamental problems in numerical magnetohydrodynamics is the formulation of 
boundary conditions. The governing equations describe the velocity and magnetic 
field in in a finite container of electrically conducting fluid. At the container bound- 
aries, the velocity is specified, but the magnetic field is not. Instead, the magnetic 
field is required to satisfy continuity conditions with the exterior magnetic field in 
the domain surrounding the fluid. The nature of these conditions depends on the 
properties of the surrounding medium; a complete discussion can be found in [il7il . 

Under the quasi-static approximation [fTTII , for a given velocity field u and magnetic 
Reynolds number Rm, the equations describing the interior magnetic field are 



5tB = V X (u X B) + -— AB (la) 

Rm 

V-B = (lb) 

The case of a fluid of finite electric conductivity restricted to a finite volume and 
surrounded by vacuum is of special importance because it models a number of ex- 
perimental, geophysical, and astrophysical configurations. Since there are no elec- 
trical currents in a vacuum, B^'"'^ is curl-free, and is therefore the gradient of a 
potential if the exterior domain is simply connected. The exterior magnetic field 
then obeys: 



B^'^^ = V0™' (2a) 
A(/)"'^^ = (2b) 
B™^-.0 |x|^cx) (2c) 

but is otherwise not fixed. The magnetic field is required to be continuous at the 
boundary: 

B - B™" = X G (3) 



In this case, continuity of all three components of the magnetic field are sufficient 
to uniquely determine both the interior and exterior fields. Our ultimate goal is 
to transform Q-© into boundary conditions that can be applied to ([T]) without 
calculating B™'^. The general principle we will employ is to construct a complete 
basis of exterior solutions B™'^ of ([21) in a preprocessing step, and to calculate B\qq 
for each member of the basis. The matching conditions Q will then yield boundary 
conditions for B. 

To explore this approach, we will apply it to the simpler analogous scalar prob- 
lem of the Poisson problem in an interior domain. We will require the solution to 
match continuously to an exterior solution satisfying Laplace's equation. In this 
case, both Dirichlet and Neumann matching conditions are necessary to specify a 
unique solution. Formally, we wish to solve the following problem: 
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A$ = p inQ (4a) 

A0 = outside Q (4b) 

with boundary conditions: 

$(x) - 0(x) = X G (91] (5a) 

a„$(x) - a„0(x) =0 X G an (5b) 

V0(x)^O |x| ^ oo (5c) 



where is a bounded domain with boundary dVl. A physical interpretation of ©- 
([5]) is that of an electrostatic potential $ of a field generated by charges distributed 
in space with the density — p, where the electrical permeability of the vacuum is 
taken to be one by the choice of units. We wish to calculate the interior solution $ 
without explicitly constructing (p. 

There exists a vast literature on the numerical solution of the fundamental physical 
problems dl])-© and ([T])-©. We will briefly survey a small portion of this literature 
here, and postpone a more detailed comparison between our method and others to 
a later section. 

The main tool by which exterior domains can be eliminated is Green's theorem, 
which replaces elliptic differential equations over a domain with integrals over the 
bounding surfaces. The use of methods based on boundary integrals has grown ex- 
plosively since the 1970s-1980s to solve engineering problems from fields such as 
acoustics, elasticity, electromagnetism and fluid mechanics HI 8 II 912012 1 122123 124125 126II . 
Hybrid methods, coupling a differential equation formulation in a domain and a 
boundary integral formulation at the boundary via an influence matrix, were also 
developed at the same time to solve (HI)-® and similar problems. The majority 
of these approaches have been based on finite elements and are hence applicable 
to complicated real- world geometries. The boundary integrals are discretized with 
techniques derived from finite element theory, leading to the term boundary ele- 
ment method, and the hybrid methods use finite elements to solve the equations in 
the domain. 

This situation contrasts with magnetohydrodynamics, which has been dominated 
by spectral methods. Spherical domains are standard, for geophysical and astro- 
physical reasons. Spectral methods can then be based on spherical harmonics and 
the poloidal-toroidal decomposition n7l8l9llOII . The solution to Laplace's equation 
on the exterior is immediate and, moreover, solutions and associated boundary con- 
ditions for each spherical harmonic and toroidal or poloidal component are decou- 
pled. Boundary conditions at the interface can then be formulated for each mode 
without the use of an influence matrix. 

The technique which we will describe is based on spectral methods, but the geome- 
try is assumed to be somewhat more complicated. Our technique occupies the niche 
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which spectral methods staked out in the 1970s- 1980s when the use of Chebyshev 
polynomials became commonplace f 27l28ll to represent domains with one or more 
non-periodic directions. In keeping with this tradition, we expect its main applica- 
tion to be to tensor-product domains whose boundaries consist of a small number 
of piecewise- smooth surfaces, such as the finite three-dimensional cylinder which 
is our eventual goal. 

We also mention here some other recent approaches to solving the magnetohy- 
drodynamic equations P 12I13I14I15I16II . with a view to generalizing the geometry 
and/or increasing parallelization. Li [I13I14II . a finite volume method is used to dis- 
cretize the solution in the interior, which is matched to that in the exterior vacuum 
via a boundary element method. (15] describes an integral equation formulation for 
the entire domain, and [|16I uses finite elements with a penalty method to apply 
boundary conditions. 



2 Influence matrix formulation 



We formulate a two-stage method for solving dH)-®, consisting of an initial pre- 
processing step which depends only on the geometry, followed by a step whereby 
solutions for many different distributions p can be generated at little incremental 
cost. This is the usual description of the decomposition of $ into homogeneous and 
particular solutions, with the additional proviso that solutions in the exterior do- 
main are to be taken into account in the preparation of the homogeneous solutions. 
We will construct the homogeneous solutions by generating harmonic bases {$^} 
and {(pj} of interior and exterior solutions, corresponding to Dirichlet boundary 
data {fj} to be specified later. We decompose (H])-® into the Poisson and Laplace 
problems: 



A$P = p inn 



^"Idn = 



A$J = in 

^^\on = fj 



A(j)j = outside Q 
(pj\dn = fj 

V0j|oo = 



(6) 



and then construct the linear superpositions: 



(j) = J2j Cj(j)j 



(7) 



Then 
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A$ = 
A(f) = 

($-0)|g^=$^|an 

V0|oo = 



j 

E^i V0j|oo = 



an 



in Q 
outside 



(8) 
(9) 
(10) 
(11) 



are already satisfied by construction, while 



an 



(12) 



constitutes a system of equations to be solved for Cj, where the derivative with 
respect to the normal is taken in the direction from the interior to the exterior region 
for both $ and 0. $ is then set equal to the sum in (|7]). If the interior harmonic 
functions are not stored, $ can be obtained by solving: 



A$ = p inn 
j 

Using Xj to index points on the boundary, (fT2l) can be discretized as: 

E[5.(0,-<f?)(x.)]c, = 9„<|.^(x,) 



(13) 
(14) 



(15) 



Equation (1151) shows that the goal of the preprocessing step is the construction and 
inversion of the influence or capacitance matrix: 



(16) 



The functions {fj} are required to constitute a complete set for values along the 
discretized boundary dCl. Another way to describe the influence matrix is as a dis- 
crete representation of the difference between the Dirichlet-to-Neumann mappings 
in the exterior and in the interior regions. 

Equivalently, 9„(0j — $^)|en can be represented as coefficients of a basis set {gi} 
(which may be identical with the set of boundary value functions {fi}) along each 
boundary. Equation (fT2l) is then discretized as: 



E(^"( 



(17) 



Although we will use and the notation in (fTSi) in what follows, the method is 
easily reformulated using (fT71) . 
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3 Solution of Poisson and Laplace problems 

3. 1 Interior domain 

We now turn to the solution of For the interior problems listed in the first two 
columns, we assume that we dispose of a solver able to compute solutions to Pois- 
son's equation in Vt with any specified boundary values. In principle, any numerical 
method can be used. In our particular case, we use a spectral discretization [i28l 

K,L 

$(x,|/)= %{x/H)%{y) (18) 

for the rectangle [—H, H] x [—1, 1]. The spectral basis functions are the Cheby- 
shev polynomials Tk{x) = cos(A; arccos(a;)). Taking > 1, we set K > L. We 
use a standard method fl^ to solve the Poisson equation with Dirichlet boundary 
conditions, diagonalizing the discretized second derivative operator in y, and using 
recursion relations to treat the second derivative in x. 

3.2 Exterior harmonic functions 

Our main focus is on the construction of the exterior harmonic solutions, specified 
in the third column of In order to avoid truncating or spatially discretizing 
the exterior domain, we will construct {(pj} using the fundamental solution of the 
Laplace equation: the Green's function satisfying 

Ax'G'(x;x') =5(x-x') (19) 
VG'(x;x')= forx'^cx) (20) 

For a specified boundary value distribution /(x), we first calculate an appropriate 
source distribution a(x) on the boundary by solving the integral equation: 

/ G'(x;xXx') = /(x) for^edn (21) 

The exterior harmonic function 0(x) required is then: 

0(x)= / G(x;x')a(x') (22) 
Jan 

where x takes values either on or off dfl. 

We now apply (|2TI)-(|22|) to our particular test problem of a rectangle. We divide the 
set of boundary distributions into four sets, each taking non-zero values on only one 



6 



side of the rectangle. In two dimensions, the fundamental Green's function solving 
(|20l)is 

^ln|x-x'| (23) 

Equation (|2T]) thus reduces to: 

/ \n\x - x'\a{x')dx' = f{x) (24) 

Ja 27r 

Equation (l24l) is known as Symm's or Carleman's equation and has the following 
solution [,29i30il : 



aix] 



{x — a)(h — x) 



b ^{x' - a){h- x')f'{x') dx' 



b f{x') dx' 



ln((6-a)/4) Ja - a)(6 - x') 



(25) 



if6 — 07^4. (If6 — a = 4 then the second integral in (|25l) can be replaced by an 
arbitrary constant if f{t) [{t - a){h - t)Y^ dt = 0.) 



Up to now, we have not specified the Dirichlet boundary values /. The choice of 
boundary value distributions is restricted only by the requirement, stated in section 
[2l that the set of distributions form a basis for functions defined on the boundary 
dfl. Because we use Chebyshev polynomials to represent the interior solutions, it is 
convenient to take as boundary values /fc(x) each of the functions T^^x/H) on the 
interval [—H, H]. The corresponding solutions ak{x; H) obtained from evaluating 
(1251) are: 



ak{x; H) = Ah 



%{x/H) 



A. 



27ik 



k>0 



-27r[ln(i7/2)]"^ A; = 



(26) 



This remarkable property - the fact that that weighted Chebyshev polynomials are 
also obtained as the source distributions ak{x) when the boundary values fk{x) 
are Chebyshev polynomials - is related to the very reason that Chebyshev poly- 
nomials are optimal in approximating polynomials on the interval. The function 
— x^ in (|26|) (for iJ = 1) is the weight with respect to which Chebyshev 
polynomials are orthogonal on the interval and is the asymptotic density of the 
Chebyshev interpolation points cos(7rj/ J), the extrema of the Chebyshev polyno- 
mials. See [3 1 32J for further details. Note also that the orthogonality of the Cheby- 
shev polynomials with respect to this weight causes the second integral in (|25l) to 
vanish except for %. 
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The corresponding harmonic functions (f)l are constructed via 



= I ^ 1^ 1^ " (^k{x'; H) dx (21) 



and are illustrated in figure [U Specifying values along along the segment [—1,1] in 
the y direction, we obtain: 

4>n^) = j\ ^ 1^ 1^ - y'^y\ ^'(^'; 1) ^y' (28) 

Notethatafc(a;;if) = ak{x/ H-l)/ H fork > and ao(a;; i/) = -^^^^ao{x/ H-l) 

The harmonic functions corresponding to specified values along the lower or upper 
boundaries {y = ±1, \x\ < /f) are 0^' (x =]= le^^); those corresponding to the left 
or right boundaries (x = rLH, \y\ < 1) are </)f (x =f Hex). We do not require the 
functions (j)%, (jy( either inside or outside Vt, but only the values and normal deriva- 
tives on the boundary. Although, for example, the values of 0^ (x + le^,) on the 
lower boundary are merely the specified values fk{x), its values on the other three 
boundaries must be calculated via (ITTI) . When evaluating the normal derivatives, 
the kernel (^(x; x') is differentiated before integration: 



-1 I 

9y(t>l{x, y) = dy — In a/ (x - x')^ + ?/2 ak{x'] H) dx' 

J-H ZTT ^ 

and similarly for S^j^f . 

Any exterior harmonic function can be approximated by the truncated series: 



K-l 



lK,L 



;x) = E [^k-<Pl (x + Hex) + Cu^<Pl (x - Hex] 



k=0 

+ E [cf'-^r (x + le,) + cf'+^r (x - le,)] (30) 

(=0 

The potential (/)(x) of (|30l) is defined by the 2(A'+L) coefficients {c^'~, c^'^, cf'~, cf' 
A very important property of the harmonic basis {(f)f.{x), 0f (x)} is that it repre- 
sents a near-boundary field uniformly. This means that the truncated series (|30l ) 
converges uniformly for any smooth boundary data and for all locations x near 
the boundary. This is a direct consequence of the excellent convergence properties 
of Chebyshev approximation applied to (|24|) : a proof can be found in [|33l . This 
property does not necessarily hold for other harmonic bases, in particular spherical 
harmonics, for which near-boundary convergence cannot be achieved, leading to a 
strong Gibbs effect. 
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In the taxonomy of boundary integral methods, equations (|2TI) and (l22l) constitute 
an indirect method, in that the intermediate surface charge density a is constructed; 
this is done by solving the Fredholm integral equation of the first kind (l24l) . The 
surface charge density is a single-layer rather than a double-layer (dipole) poten- 
tial; equivalently G, rather than dG/dn, is used in the representation. Because the 
method determines only the boundary values and normal derivatives of the exte- 
rior solution, it is not the preferred approach when the exterior potential is itself 
required at each time step: although the exterior solution can be sampled at any lo- 
cation, this is computationally expensive, as is often the case for boundary integral 
methods. 




Fig. 1. Potentials (/)^(x) generated by line source distributions ak{x;H = 2). Values of 
the potentials on the line segment x G [—1,1] correspond to the Chebyshev polynomials 
Tk{x). From top to bottom: (To(x), o"i(x), ai{x). 
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4 Electrostatic example 



We apply our method to a simple problem from electrostatics, the distribution of 
electric charges ~Pm confined in a rectangular domain but localized around the 
origin: 



r^e ^'^/''^ cos{m9) for |x| < 1 and \y\ < 1 







for > 1 or |y| > 1 



(31) 



The factor in (|3T|) ensures regularity of at r = 0. The potential due to 
unbounded sources (not restricted to the interior domain) can be found analytically: 



$m=i(r,( 



4 

4r 
4r2 



Edl,- +21og(r) 



e ^ 



2 



cos 



{5' + 



) e 



r2 



COS 2^ 



(32a) 

(32b) 
(32c) 



where Ei{a, A) is the error function Ei{a, A) = j"^ e ^^r " rfr and are chosen 
to be finite at r = 0. 

We seek the corresponding electric potential. The parameter S is chosen to make pm 
very small near the boundaries. We expect the solution to be almost unaffected by 
the presence of boundaries. The source distribution Pm=o (r, 6) should therefore lead 
to a potential which is almost axisymmetric. Figure [2] shows the potential obtained 
numerically for 5^ = 0.15 using the spectral resolution = 8 in both directions. 
The domain boundary is represented by a bold square. The contours are almost 
perfectly circular, as should be the case for 5 small, showing that the presence of 
the boundaries has minimal effect. 

To evaluate the error convergence of the method we computed the relative error 
E„i{N) defined as 

\^Ur,e)-^Zir,e)\ 

E„\N) = sup 1^ . (33) 

where $^(r, 6^) is the solution computed numerically with spectral resolution N 
in both spatial directions of the bounding square and ^mi^, 6) is the analytic so- 
lution (l32l) in the absence of the bounding square. Figure [3] proves the exponential 
convergence of the method. 

Figures |4]-[5] show the electric potentials and $^=2^ for 5 = 0.1. Figure 

[6] shows $m=i^ with (5 = 2. In figure Ul the dipole source distribution has been 
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Fig. 2. Potential ^>^=^ for 
6 = 0.15. Maximal relative error is 



Fig. 3. Convergence test: Conver- 
gence test: log^Q £'m,=o(^) plotted for 
N = [6,...,16],(52 = 0.1. 



Em=o{N = 8) « 0.03. 



rotated by 45° about the origin. For this large value of 6, charges are located near 
the boundary. In each case with 5 <^ 1, we observed exponential convergence 
toward solution (|32|) . Convergence can only be confirmed up to a limited precision 
since the analytic solution (|32l) does not correspond exactly to the problem we are 
solving numerically, in which sources are confined to the interior square. The best 
agreement can be achieved for small values of S. If the numerical solution with 
highest spectral resolution (here = 64) is instead taken as a reference, then the 
method converges to this solution spectrally up to machine precision for any value 
of 5. 



5 Implementation 

5. 1 Summary and computation cost 

We describe the implementation of the method for our illustrative example (fTSi) of 
the rectangle [—H, H] x [—1,1] with double Chebyshev discretization (i^' + 1) x 



The total preprocessing step consists of: 

• Evaluation of the values and the normal derivatives of the exterior harmonic so- 
lutions on the boundary. 

• Calculation of the interior harmonic solutions. 

• Inversion or LU decomposition of the influence matrix, 

For each particular right-hand-side p, the operations consist of: 
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Fig. 6. Potential ^^=1^ for S = 2. Fig. 7. <^^zf{r, + f ) for 5 = 2. 

• Solving for a single particular solution. 

• Acting with the inverse of the influence matrix. 

• Using the corrected Dirichlet boundary conditions to calculate the final solution. 

The total number J of boundary points is 2{K+L). The inversion orLU-decomposition 
of the influence matrix C in the preprocessing stage requires a time proportional to 
J^, while each solution of the linear system (fT5l) determining the coefficients of the 
homogeneous solutions requires a time proportional to J^. Each interior solution is 
calculated at a cost proportional to KL'^. 

Symmetry can be used to reduce the cost of each step. The symmetry of the rect- 
angle divides all the independent harmonic solutions into four mutually orthogonal 
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Calculation 


Result Cost 


Preprocessing 


Exterior harmonic solutions 

Interior harmonic solutions 

Influence matrix inversion/decomposition 


Kx (Int, + Int^) 

(Pj\dQ,, On<Pj\dn 

+L X (Inty + Int^) 
dn<^^ (K + L)KLy4: 
C7-1 {K + Lf/2 


For each right-hand-side 


Particular solution 

Action with influence matrix 

Corrected solution 


KL'^/2 
Cj {K + Lf 
$ KL'^/2 



Table 1 



Operation count of each preprocessing and right-hand-side-dependent step for a rectangle 
discretized with {K + 1) x (L + 1) Chebyshev polynomials and points. Int^^ and Int^ are 
the costs of performing the singular integrals over x in (l27l) and ( |29l ). and Inty and Int^ are 
those of the analogous integrals over y. 

classes. Decoupling the Laplacian operator according to parity in x and y leads to 
four Poisson problems, each with resolution K/2 x L/2, thus reducing the time by 
a factor of two. Decoupling by parity also reduces the influence matrix C to four 
matrices, the dimensions of which are one fourth of that of the original matrix. 

Table 15. II gives the operation count of each step, taking into account the reductions 
permitted by symmetry. 

5.2 Singular integrals 

The integrations in (l27l)- (|29l) are performed numerically. Special attention must be 
paid in doing so since both the kernel G'(x; x') and the density cr(x') have integrable 
singularities within the domain of integration. The singular points are x'e^ = x and 
x' = ±H for (pl and y'ey = x and y' = ±1 for 0f . Dedicated adaptive quadratures 
(see Il34ll ) can be used to compute these integrals accurately. 

It is also possible to evaluate the singular part of the integral analytically, reduc- 
ing the numerical problem to the evaluation of integrals with non-singular inte- 
grands. The remaining integrand is piecewise C°° and can be integrated with spec- 
tral precision over each of the regular subdomains. Singularity subtraction greatly 
decreases the variation in grid density needed to sample the integrand homoge- 
neously, thereby significantly accelerating the numerical integrations in (l27l)-(|29l). 
Specifically, an adaptive method requires a smaller number of iterations, or, alter- 
natively, a non-adaptive method requires a coarser resolution. However, the con- 
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Fig. 8. Quadratic fit of the condition number C{N) of the influence matrix defined in ( fT6l ). 

vergence of our approximation with K, L is exponential (spectral), regardless of 
whether the singular part of the integral is subtracted or included in the numerical 
evaluation. 



5.3 Conditioning of matrices 

The influence matrix (fT6l ) is not immediately invertible. Because of the redundancy 
of information at the comers, this matrix has exactly four zero singular values or 
eigenvalues. The corresponding linear system can be solved after arbitrarily cor- 
recting singular values or eigenvalues of the influence matrix; see f,35.36J for more 
details. 

The condition number C of the corrected matrix depends on the spatial resolution 
and the maximal order of derivatives used to express the boundary conditions. In 
our case of a Neumann boundary condition and a resolution N in each direction, the 
condition number scales as C = O (iV^). Fitting the condition numbers computed 
for G [2, 32] with a parabola (see fig. 15.31) yields a formula for predicting the 
condition number for an arbitrary resolution: 

C{N) = 3.58N^ - 18. 16N + 59 (34) 

It can then be deduced from (|34l) that a reasonably conditioned matrix with C < 10^ 
is obtained for a spatial resolution as high as 0(1000). The method can therefore 
be applied to problems where small-scale field features require use of high spatial 
resolution. 



C(N) 




3.58iV2 _ 18.167V + 59 
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6 Generalizations 



We now discuss the applicability of this method to other geometries, problems, 
and spatial discretizations. The decomposition into interior particular and homoge- 
neous functions and exterior homogeneous functions described in section [2] is, of 
course, completely general and not related to any particular spatial discretization. 
The method described in section 13.11 for constructing the exterior harmonic func- 
tions relies on the Chebyshev-Chebyshev discretization of the rectangle which is 
widely used since the Chebyshev polynomials are optimal approximants of smooth 
functions. This property, as well as the straightfoward correspondence between in- 
terior and exterior solutions, make the Chebyshev discretization especially suitable 
for the construction of the exterior harmonic solutions as well. However, the method 
is easily generalizable to other basis functions /^(x) for the potential values, which 
can be substituted into (|25l) in order to calculate the corresponding charge densi- 
ties crfc(x), if dictated by the geometry or numerical method used for the interior 
problem. 

In fact, since our real interest is in generating the complete set of crfc(x) necessary 
to generate the complete set of (j)k{x), rather than in calculating the specific ak{x) 
corresponding to each particular fk{x), the only information really required in (|26|) 
is the singularity l/vrVl — x'^. One may then allow the set of cr^'s to be the prod- 
ucts of this singularity with the members of any appropriate basis set of analytic 
functions on the boundary in question. 

In three dimensions, the fundamental Green's function is 

G(x,x') = -^^- (35) 
47r |x — x'l 

In an axisymmetric geometry with a Fourier representation of the azimuthal direc- 
tion, all of the problems to be solved decouple according to Fourier mode. The 
operation count would then scale linearly with the number of azimuthal points or 
Fourier modes. The elliptic problems in Q would remain two-dimensional, and the 
integral equations equivalent to (l24l) would remain one-dimensional. 

This method can also be applied to other elliptic problems or to parabolic problems. 
As stated in section [H our motivation for developing this method is to apply it to 
the magnetohydrodynamic equations Q-©, in which ([T]) is a parabolic equation. 
A general parabolic problem may be written as: 

dt<^ = m + J^{^,p) (36) 

where T may include nonlinear and/or time-dependent source terms. First-order 
implicit temporal discretization of (l36l) results in the inhomogeneous Helmholtz 
equation: 

(/ - 5tl^)^{t + 5t)=J^ (37) 
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where T may depend on previous values of $. This Helmholtz operator (J — 
(5tA) can replace the Laplacian in (01) and (|6l). It is known that replacement of 
the Helmholtz equation by a boundary integral equation can lead to singulari- 
ties for certain values of the wavenumber (here i/\f5i); a large body of work, 
e.g. B7I38I39I40I411I . addresses this problem. However, in the magnetohydrody- 
namic case of a conducting fluid surrounded by an exterior vacuum, no such dif- 
ficulties would be introduced, since the exterior problem remains governed by 
Laplace's equation. More complicated vectorial operators may appear, as occur in 
the Navier-Stokes or magnetohydrodynamic equations. 



7 Comparison with other approaches 

We mention here some other techniques that have been used to solve exterior prob- 
lems or to match interior and exterior domains. Spectral methods can be combined 
with various transformations and mappings. The inner region can be surrounded by 
a sphere, and the outer domain decomposed into the region inside and outside the 
sphere. The exterior domain can be mapped into an interior domain via a 1/r map- 
ping [16142143 1: spectral methods can then be used to treat either or both domains. 
The region exterior to one or more spheres has been mapped to the interior of a 
rectangle [|44l or a pentangular [|45l region rotated about an axis, and Chebyshev- 
Fourier expansions used to solve elliptic equations arising in the study of black 
holes in general relativity. A smooth boundary can be parameterized by angle, and 
the boundary values represented as a series of trigonometric functions or spherical 
harmonics 112213914014 1142143146147148)1 . Our method differs from these in that a 
Chebyshev approximation is used to represent the boundary values on each seg- 
ment of a non-smooth boundary, and an analytic formula is used to calculate the 
surface density which exactly yields this Chebyshev approximation. 

Conformal mapping is another technique which can be used to calculate interior or 
exterior harmonic functions. The Riemann mapping theorem guarantees the exis- 
tence of a conformal transformation from the interior or exterior of a simply con- 
nected domain into the interior or exterior of a unit disk; its proof is, however, 
non-constructive, and does not explicitly derive the transformation. For some ge- 
ometries, including the exterior of a rectangle, an analytical formula can be derived 
[|49l . For polygon-bounded regions with piecewise-constant boundary conditions, 
the Schwarz-Christoffel [50| mapping has proved to be a very robust tool, applied 
to problems arising in magneto- and electro- statics, potential flows, inverse prob- 
lems and many other fields. 

Our influence matrix approach relies on calculating harmonic functions with ar- 
bitrary boundary data, for which conformal mapping is much more problematic. 
More general conformal mappings are often computed by solving Symm's or Car- 
leman's equation (|24l) numerically on the domain boundary, making this approach 
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similar in terms of numerical cost and precision to boundary integral equation meth- 
ods. It is interesting to note that, for domains including corners, the Chebyshev ap- 
proximation is especially well suited, guaranteeing superconvergence of the map- 
ping function [|33l . 



8 Conclusion 

As a test case for the magnetohydrodynamic equations, we have developed a method 
for solving the two-dimensional Poisson equation in a bounded domain, where the 
solution satisfies matching conditions with a harmonic potential outside the do- 
main. The method solves only the interior problem and determines the boundary 
conditions ensuring smooth matching with the exterior solution. The essential ele- 
ment of this approach is construction of a basis of harmonic functions which repre- 
sent the near-boundary exterior solutions uniformly. This basis is used to construct 
the influence matrix which serves to impose the coupled boundary conditions be- 
tween the interior and exterior solutions. The method is numerically reasonably 
well conditioned and can be used for quite high spatial resolutions. For a spectral 
solver, this method guarantees exponential convergence. 

Instead of corresponding to point sources on the boundary, each exterior harmonic 
solution corresponds to a spectral basis function. The most costly process - the 
construction of a basis of exterior harmonic functions - depends only on geometry 
and spatial resolution. Once the basis is computed it can be stored and used for any 
computation using the same resolution and domain shape. When used as a prepro- 
cessing step for time-dependent simulations, the cost of constructing the exterior 
harmonic basis is negligible compared to that of thousands of time steps. 

Since cylindrical coordinates have one periodic direction, it should be possible to 
apply this method separately to each of the Fourier modes, treated individually as 
two-dimensional problems. The extension of this method to the MHD equations in 
a finite cylindrical geometry is currently under investigation. 
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